library(tidyverse)
library(lfe)
N <- 100

# Time period data
tp1 <- data.frame(exposure = rnorm(N),
                  error = rnorm(N),
                  unit = 1:N,
                  time = 1)
tp2 <- data.frame(exposure = rnorm(N),
                  error = rnorm(N),
                  unit = 1:N,
                  time = 2)
tp3 <- data.frame(exposure = rnorm(N),
                  error = rnorm(N),
                  unit = 1:N,
                  time = 3)

df <- bind_rows(tp1, tp2, tp3)

# Add in unit effects and time effects
unit_effects <- data.frame(unit = 1:N, unit_effect = rnorm(N))
time_effects <- data.frame(time = 1:3, time_effect = rnorm(3))

# Craate large df
df <- df %>%
  left_join(unit_effects) %>%
  left_join(time_effects)

# Create outcome
df <- df %>%
  mutate(
    beta = case_when(time == 1~1, time == 2~3, time == 3~8),
    y = unit_effect + time_effect + error + exposure*beta
  )

# Try regressions
spec1 <- felm(y ~ I((time == 2)*exposure) + I((time == 3)*exposure)  | 
                unit_effect + time_effect | 0 | unit_effect,
              data = df)
spec2 <- felm(y ~ exposure +  I((time == 2)*exposure) + I((time == 3)*exposure)  | 
                unit_effect + time_effect | 0 | unit_effect,
              data = df)
spec3 <- felm(y ~ I((time == 1)*exposure) +  I((time == 2)*exposure) + I((time == 3)*exposure)  | 
                unit_effect + time_effect | 0 | unit_effect,
              data = df)